#!/usr/bin/python
import Gnuplot
import math
from random import uniform


def markov_throw(delta, x, y):
	return	x + uniform(-delta, delta), y + uniform(-delta, delta) 


def direct_throw(start, end):
	return uniform(start, end), uniform(start, end)


def IsInSquare(x, y):
	if abs(x) < 1. and abs(y) < 1.:
		return True 
	else:
		return False
	

def IsInCircle(x, y):
	if x**2 + y**2 < 1.:
		return True
	else:
		return False


def CalErr(summ, sqrSum, N):
	avg = summ/float(N)
	st_err = (sqrSum/float(N)-avg**2)**1/2
	return st_err


def markov_pi(delta, N):
	x, y = 1., 1. # initial position
	N_hits = 0
	record = []
	sumPi = 0.0	#sum(f)
	sumSqrPi = 0.0	#sum(f^2)

	for i in range(1,N):
		tx, ty = markov_throw(delta, x, y)

		if IsInSquare(tx, ty):
			x, y = tx, ty

		if IsInCircle(x, y):
			N_hits += 1

		c_pi = N_hits * 4. / float(i)
		sumPi += c_pi
		sumSqrPi += c_pi**2

		if i > 2:
			err = CalErr(sumPi, sumSqrPi, i)
		else:
			err = 0.0

		if i % 100 == 1:
			record.append([i, c_pi, err])

	return record


if __name__== "__main__":
	#data = markov_pi(.3, 10001)
	#item_pi = Gnuplot.PlotItems.Data(data, cols=(0,1), title = 'pi', with_='lp pt 6')
	#item_err = Gnuplot.PlotItems.Data(data, cols=(0,2), title = 'err', with_='lp pt 6')
	#g1 = Gnuplot.Gnuplot()
	#g1('set multiplot title "Markov_pi" layout 2,1')
	#g1.plot(item_pi)
	#g1.plot(item_err)

	N = 10001 # iter times
	delta_value = [0.01, 0.1, 0.3, 0.5, 0.99] 
	pi_items = []
	err_items = []

	for i in delta_value:
		calc = markov_pi(i, N)
		pi_items.append(Gnuplot.PlotItems.Data(calc, cols=(0, 1), title='pi '+str(i)))
		err_items.append(Gnuplot.PlotItems.Data(calc, cols=(0, 2), title='err '+str(i)))

	g_pi = Gnuplot.Gnuplot()
	g_pi.xlabel('N')
	g_pi.ylabel('pi')
	g_pi('set yrange [0:4.7]')

	g_err = Gnuplot.Gnuplot()
	g_err('set logscale xy')
	g_err.xlabel('N')
	g_err.ylabel('err')

	for i in range(len(delta_value)):
		g_pi.replot(pi_items[i])
		g_err.replot(err_items[i])

	raw_input('Press Enter to exit')
